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The effective interaction between charged colloidal particles confined between two planar like- 
charged walls is investigated using computer simulations of the primitive model describing asym- 
metric electrolytes. In detail, we calculate the effective force acting onto a single macroion and onto 
a macroion pair in the presence of slit-like confinement. For moderate Coulomb coupling, we find 
that this force is repulsive. Under strong coupling conditions, however, the sign of the force depends 
on the distance to the plates and on the interparticle distance. In particular, the particle-plate 
interaction becomes strongly attractive for small distances which may explain the occurrence of 
colloidal crystalline layers near the plates observed in recent experiments. 

PACS: 82.70.Dd, 61 .20.Ja 



I. INTRODUCTION 

There is recent experimental evidence that the effec- 
tive interaction between like-charged colloidal particles 
( "macroions" ) is sensitively affected by a confinement be- 
tween two parallel charged glass plates For aque- 
ous polystyrene suspensions studied in experiment, the 
effective force between two colloidal macroions is found 
to be repulsive far away from the plates but becomes at- 
tractive when the like-charge macroions are located close 
to an equally charged plate. At first glance, these find- 
ings are surprising as one would expect a purely repulsive 
interaction from the electrostatic part of the traditional 
Derjaguin-Landau-Verwey-Overbeek (DLVO) theory |Q. 
In fact a full theoretical explanation is still missing but 
several steps were performed in different directions: the 
essential difference in a confining geometry with respect 
to the bulk is that the counterion density field is inho- 
mogeneous for small coupling between the macroions and 
counterions. In a straightforward generalization of the 
DLVO theory to such an inhomogeneous situation H|), 
the effective force between the macroions remains repul- 
sive close to the charged plates but becomes weaker since 
the local concentration of counterions is higher which re- 
sults in a stronger screening of the Coulomb repulsion. It 
was further realized that a charged wall induces signifi- 
cant effective triplet interactions |7]] which are ignored in 
the usual DLVO approach resulting in a net attraction 
Q or in a repulsion |^] depending on the system param- 
eters. An explicit calculation was done within density 
functional perturbation theory which is justified, how- 
ever, only for weak inhomogeneities. A complementary 
approach is to solve the nonlinear Poisson-Boltzmann 
equation with appropriate boundary conditions in a fi- 
nite geometry. This was done recently for two charged 
spheres in a charged cylindrical pore J10| as well as for 



two charged cylinders confined by two parallel charged 
plates |11|. However, the first situation corresponds to 
a finite system where the Poisson-Boltzmann approach 
does not lead to attraction Jl^] and the second situa- 
tion is a quasi-two-dimensional set-up which is known to 
behave qualitatively different from a three-dimensional 
situation [ jl3| . A further complication arises from image 
charges induced by the different dielectric constants of 
the glass and the solvent Jl^l6[|. 

A general problem of any theoretical description (as 
DLVO, Poisson-Boltzmann) is that close to the walls the 
counterion concentration is high and any weak-coupling 
theory fails a priori when applied to a situation of con- 
fined macroions. For strong coupling, even in the bulk, it 
is unclear whether an effective attraction of like-charged 
spherical macroions is possible although there are hints 
from experiment s jl^ |l9|], theory [p0[-p4|, and computer 
simulations [p5^p7fl . At this stage it is important to re- 
mark that a phase separation seen in experiment does 
not necessarily imply an effective attraction. The ad- 
ditional contribution from the counterions to the total 
free energy may induce such a phase separation although 
the effective interaction between the macroions is purely 
repulsive ]28| , p9| ]. Bearing the difficulties in experimen- 
tal interpretations and theory in mind, computer sim- 
ulations represent a helpful alternative tool to extract 
"exact" results for certain model systems. The general 
accepted theoretical model for the description of charged 
colloidal suspensions is the "primitive approach" where 
the discrete structure of the solvent is disregarded and 
the interaction between the macroions and counterions is 
modelled by excluded volume and Coulomb forces. The 
problem with a full computer simulation of the primitive 
model is the high charge asymmetry between macro- and 
counterions which restricts the full treatment to micelles 
rather than charged colloidal suspensions [flof . 

In this paper, we use computer simulations to obtain 
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"exact" results for the effective interaction between con- 
fined charged colloids based upon the primitive model. 
Instead of solving the full many-body problem with many 
macroions, we only simulate one or two macroions con- 
fined between two parallel charged plates. This enables 
us to access high charge numbers of the macro-particles. 
As a result we find that the wall-particle and the inter- 
particle interaction is repulsive for weak Coulomb cou- 
pling. For stronger coupling, the behaviour of the force 
changes from repulsive to attractive and back to repul- 
sive as the interparticle distance is varied. In particular, 
the plate-particle interaction exhibits a short-range at- 
traction for a small distances. This may explain the oc- 
currence of crystalline colloidal layers on top of the glass 
plates found in recent experiments |3l|-j33l|. These crys- 
tallites are metastable but very long-lived and cannot be 
understood in terms of DLVO-theory. 

The paper is organized as follows: the model and our 
target quantities are defined in section II. Section III 
contains details of our computer simulation procedure. 
Results for the counterion density profiles are shown in 
section IV. The case of a single macroion is discussed in 
section V, and a macroion pair is investigated in section 
VI. Finally, we conclude in section VII. 



II. THE MODEL AND TARGET QUANTITIES 

We consider N m macroions with bare charge q rn 
'/.< > (e > denoting the elementary charge) and 
mesoscopic diameter d m confined between two parallel 
plates that carry surface charge densities <j\ and 02- 
We assume that the plates and the macroions are likely 
charged. The separation distance between plates is 2L. 
For convenience, we choose the z axis to be perpendic- 
ular to the plate surface. The origin of the coordinate 
system is located on the surface of one plate. Image 
charges are neglected, i.e. we assumed for simplicity that 
the dielectric constants of the solvent, the plate and the 
colloidal material are the same. Typically we use a peri- 
odically repeated square cell in x and y direction which 
possesses an area S p . Hence the macroion number den- 
sity is p m — N m /2LS p . We restrict our studies to a small 
number of macroions in the cell. In particular we are con- 
sidering the cases N m = 0,1,2 subsequently. Both the 
macroions and the charged plates provide neutralizing 
counterions which are dissolved in a solvent of dielectric 
constant e. The counterions have a microscopic diameter 
d c and carry an opposite charge q c — ~qe where q > de- 
notes the valency. Typically, q = 1,2. For simplicity, we 
assume that the counterions from the walls and from the 
macroions are not distinguishable. The total counterion 
number N c in the cell (as well as the averaged counterion 
number density p c — N C /2LS P ) is fixed by the condition 
of global charge neutrality 
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The interactions between the particles are described 
within the framework of the primitive model. We assume 
the following pair interaction potentials V mm (r), V mc (r), 
V cc (r) between macroions and counterions, r denoting 
the corresponding interparticle distance: 



V mm (r) = < z 



00 

V mc (r) = \ Zqe 



OO 
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for r < d n 
for r > d n 



for r < {d m + d c )/2 
for r > (d m + d c )/2 



V cc (r) 



00 for r < d c 
= i £lL for r > d 
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The interaction between the particles and the wall is de- 
scribed by the potential energy 

!oo for z < di/2 and 

z > 2L - dj2 (5) 
2 " (CT2 7 l)g ' z else 

where z is the altitude of the particle center and i — 
m, c. Note that the interaction between the wall and the 
particles is zero for equally charged plates. 

Our target quantities are the equilibrium counterion 
profiles and the effective forces exerted on the macroions. 
The counterionic density profile p^ (r) is defined as sta- 
tistical average via 



(G) 



where {fj = {xj,yj,Zj);j = 1,...N C } denote the coun- 
terion positions. The canonical average < ... > c over 
an {fj }-dependent quantity A is defined via the classical 
trace 



x„4({r fc })exp 
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where ksT is the thermal energy (fee denoting Boltz- 
mann's constant) and 

^ = ^^ui-^ I) 

n=l j=l 

N c N c 

+ 9 E Vccdn-rj |) + E^fe) ( g ) 
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is the total counterionic part of the potential energy pro- 
vided the macroions are at positions 
{Rj = (Xj, Yj, Zj);j — 1, ...N m }. Furthermore, the clas- 
sical partition function 
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Z=— d 3 n... / d 3 r Nc exp 

J"c! JV JV 



Vc 



(9) 



guarantees the correct normalization < 1 > c = 1. Note 
that the counterionic density profile (r) depends para- 
metrically on the macroion positions {-Rj}. 

The total effective force Fj acting onto the jth 
macroion contains three different parts 34 35|p5[ | 



F — F y 
(i) 
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The first term, Fj , is the direct Coulomb repulsion 
stemming from neighboring macroions and the plates 




V mm {\ Ri - Rj |) + Vpm^Zj) (11) 



The second part Fj involves the electric part of the 
counterion-macroion interaction and has the statistical 
definition 
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Finally, the third term Fj describes a depletion (or con- 
tact) force arising from the hard-sphere part in V mc (r), 
which can be expressed as an integral over the surface Sj 
of the jth macroion 

i? (3) =k B T J^df pW(f) (13) 

where / is a surface vector pointing towards the macroion 
center. This depiction term is usually neglected in any 
DLVO or Poisson-Boltzmann treatment but becomes ac- 
tually important for strong macroion-counterion cou- 
pling. We define the strength of Coulomb coupling via 
the dimensionless coupling parameter [p5[ 



r — 

-L m.c. — 



Z 2\ B 

q d m + d c ' 



(14) 



where the Bjerrum length is — q 2 e 2 /eksT. 

A further interesting quantity is the countcrion- 
averaged total potential energy defined as 



U({Rj}) = 



N m 

£ 

i,j;i<j 



V mm {\ Ri-Ri\)+< V e > e (15) 



In general the effective force (|T(]) is different from the 
gradient of U({Rj}) M i.e. 



III. DETAILS OF THE COMPUTER 
SIMULATION 

The Coulomb interactions involved in the primitive 
model are long-ranged but the periodically repeated sys- 
tem is finite which poses a computational problem. This 
can be solved in different ways. The simplest way to solve 
the problem is to cut off the range of the Coulomb inter- 
action by half of the system size which is the minimum 
image convention (MIC). The MIC is easy to implement 
but serious cut-off errors can be introduced. A better 
way is to include JV periodically repeated images (PRI) 
of neighbour cells in x and y direction. Also the limit 
J\f — > oo can be treated by a suitable generalization of the 
traditional Ewald summation technique 37 39] to a two- 
dimensional system. A straightforward generalization, 
however, leads to quite massive computational effort p0[ . 
A much more effective alternative is the so-called Lekner 
summation method |4l| , fi"2"f which has recently been ap- 
plied successfully to the problem of effective interactions 
between rodlike polyelectrolytes and like-charged planar 
surfaces p3| . 

A completely different way out of the problem is to 
study the system on a surface of a four-dimensional (4D) 
hypersphere which itself is a compact closed geometry 
with spherical boundary conditions p4| . Then one has 
to express the Coulomb forces in terms of the appropri- 
ate 4D spherical coordinates which can be done analyt- 
ically, see Appendix A. Such spherical boundary condi- 
tions were effectively utilized in computer simulations of 
two-dimensional (2D) classical electrons |^,^6| and other 
2D fluids H) . Simulations of the 3D system located 
on the surface of a 4D hypersphere were carried out for 
Lennard- Jones |50) , hard sphere J5l| and charged |52j sys- 
tems. The hypersphere geometry (HSG) was also tested 
against Ewald summations to investigate the stability of 
charged interfaces ]53| and good agreement was found, 
even for strongly coupled interfaces. Simulations in HSG 
are much faster than that for Lekner sums or PRI as 
there is no sum over images. 

In most of our investigations we have used HSG simula- 
tions but tested them against MIC, PRI and Lekner sum- 
mations. Good agreement was found except for the MIC 
which suffers from the early truncation of the Coulomb 
tail. We have performed Molecular Dynamics (MD) sim- 
ulations at room temperature T — 293° K. A more de- 
tailed description of the MD procedure in HSG is given 
in Appendix B. The width of planar slit is fixed to be 
2L = 5d m . Different sets of system parameters are sum- 
marized in Table I. 



(16) 



In fact, as we shall show below these two quantities be- 
have qualitatively different for strong coupling. We em- 
phasize that it is the effective force (|l^) that is probed 
in experiments. 
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TABLE I. Set of parameters used in our calculations 



Run 


N m 


Z 


q 


<7 p (e/cm 2 ) 


e 


Pm(l/cm 3 ) 


d m (crri) 


d c (cm) 




A 





- 


2 


1.24 x 10 11 


78.3 


1.17 x 10 i3 


- 


5.32 x 10" 8 


- 


B 


1 


200 


2 


0.62 x 10 11 


78.3 


1.17 x 10 13 


5.32 x 10" 6 


5.32 x 10" 8 


11 


C 


1 


200 


2 


1.24 x 10 11 


78.3 


varied 


5.32 x 10~ 6 


5.32 x 10~ 8 


11 


jj 


I 




2 


1.49 x 10 11 


1 1/1 f"! fl 


1.17 x 10 13 


O .O^ A ± U 


O.iJj^ A 1U 


1 1/1 T' 'J O /7 

(JLLI fcoLl 


E 


1 


100 


2 


2.98 x 10 11 


3.9 


1.17 x 10 13 


5.32 x 10~ 6 


5.32 x 10~ 8 


110 


G 


1 


100 


2 


varied 


78.3 


9.36 x 10 1B 


2.66 x 10" 7 


2.66 x 10" 8 


100 


K 


1 


32 


2 


1.56 x 10 14 


77.3 


1.9 x 10 18 


2.56 x 10" 7 


2.56 x 10" 9 


37 


L 


2 


200 


2 


1.24 x 10 11 


78.3 


2.34 x 10 13 


5.32 x 10 -6 


5.32 x 10 -8 


11 


M 


2 


100 


2 


varied 


3.9 


2.34 x 10 13 


5.32 x 10~ 6 


5.32 x 10~ 8 


110 


N 


2 


100 


2 


varied 


78.3 


1.87 x 10 17 


2.66 x 10~ 7 


2.66 x 10~ 8 


100 



We take divalent counterions throughout our investiga- 
tions. The dielectric constant is that for water at room 
temperature (e = 78.3) but we have also investigated 
cases where e is smaller in order to enhance the Coulomb 
coupling formally. The charge asymmetry Z/q ranges 
from 16 to 100. The time step At of the simulation was 
typically chosen to be 10~ 3 \J m d^/e 2 , (with m denot- 
ing the mass of the counterions) such that the reflection 
of counterions following the collision with the surface of 
macroions and walls is calculated with high precision. For 
every run the equilibrium state of the system was checked 
during the simulation time. This was done by monitoring 
the temperature, average velocity and the distribution 
function of velocities and total potential energy of the 
system. On average it took about 10 4 MD steps to get 
into equilibrium. Then during 5 • 10 4 — 5 • 10 5 time steps, 
we gathered statistics to perform the canonical averages 
for calculated quantities. 



IV. COUNTERION DENSITY PROFILES 
BETWEEN CHARGED PLATES 

First, as a reference case, let us discuss the situation 
without any macroion. This set-up is well-studied in the 
literature |54|,|55|]. We consider equally charged surfaces 
(j i = (72 = o v . The imbalance in the interaction with 
neighbours will push the counterions toward the plates. 
Consequently, a great majority of neutralizing counteri- 
ons reside within a thin surface layer. For strong cou- 
pling, the width of this layer can be approximately ex- 
pressed as M 



\2 

z 2L' 

where A_d is the bulk Debye screening length 

ek B T 



\ 2 - 
A D — 



2„2 



(17) 



(18) 



where po = p c . Due to symmetry, the equilibrium coun- 
terion density profile only depends on z. The analytical 



solution of the Poisson-Boltzmann (PB) equation for this 
profile is |57l] 



pi PB) {z)=P, 



LW( 7o (W)) 



(19) 



where 70 is defined via the solution of the implicit equa- 
tion 



{L/\d? 
270 



tan 70 = 



(20) 



For parameters of moderate Coulomb coupling (run A), 
the PB result is shown as a solid line in Fig.[l| 
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FIG. 1. Reduced counterion density profiles p c (z)/po ver- 
sus reduced distance z/L. solid line- PB prediction and sim- 
ulation result with incorporating Lekner summation method. 
Both data do coincide on the same curve, long-dashed line- 
result of simulation in HSG, dashed line- result of PRI simu- 
lation with M = 2, dotted line- result of MIC simulation. 

The corresponding MD simulation data were obtained 
with 600 counterions in the cell using various bound- 
ary conditions. As expected, the PB theory coincides 
with the Lekner summation method which treats best 
the long-range nature of the Coulomb interactions. In 
HSG the counterionic profiles are also very similar to the 
Lekner summation while the MIC deviates significantly. 
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The MIC can already be improved significantly if TV = 2 
periodic repeated images are included. In conclusion, the 
agreement between Lekner summation and HSG justifies 
the HSG a posteriori and gives evidence that the HSG 
produces reliable results also for stronger couplings. 



V. SINGLE MACROION BETWEEN CHARGED 
PLATES 

Let us now consider a single macroion in the inter- 
lamellar area. We put the macroion on the z-axis, such 
that its position is at R\ — (0,0, Zi). A corresponding 
schematic picture is given in Fig.|^. 



FIG. 2. Schematic picture for a single macroion between 
likely charged planes of charge density a p , separated by dis- 
tance 2L. 



o 







1 



0.5 

Zj/L 

FIG. 3. Force FY = FY ■ e z acting on a single macro-ion 
versus reduced macro-ion distance Zi/L. The force is scaled 
by the (arbitrary) unit Fo = -^p— • The system parameters 
are from run B. The solid line is the prediction from PB 
theory. The open circles are simulation results in HSG. The 
statistical error corresponds to the symbol size. 

This justifies the theoretical conclusions drawn in Refs. 
!||| based on PB theory. The deviation between the- 
ory and simulation are larger in Fig.^ where the surface 
charge density was doubled. Here, also the system size 
dependence (resp. the dependence on the macroion den- 
sity) was studied in the simulation. 



The total force acting on the macroion only depends 
on Z\ and points along the z-axis. Obviously, for the 
case <7i = (72 = <J P of symmetric plates considered here, 

the direct part of the total force, Fj- , vanishes. For the 
second (electrostatic) part, simple PB-theory applied to 
the case of small macroion charge and small macroion 
diameter yields the followinganalytical expression for the 



effective macroion force 57 



PB 



2Zk B T l0 
qL 



tan (70 (1 - z/L))e z 



(21) 



where e z is the unit vector in z-dircction and 70 is given 
by ( |2p| ) . Note that only the counterion density stemming 
from the charged plates has to be inserted in ( pc| ) i.e. 

This force pushes the macroion towards the 



Po 



L\qe\ 



mid-plane, i.e. the wall-particle interaction is repulsive. 

The expression (21) will break down, however, for a 
large macroion diameter d m and for strong macroion- 
countcrion coupling parameter T mc . We have tested the 
PB-theory against "exact" computer simulation data. 
For moderate couplings (run B and run C) the results 
for the total force F\ = F% ■ e z are shown in Figs^-Q. In 
Fig.^, a surprising agreement between theory and simula- 
tion is obtained despite the fact that the macroion charge 
is large. 




FIG. 4. Same as FigJ^ but now for run C. Symbols are sim- 
ulation results in HSG for various macroion number densities: 



circles: 



pm 



1.17xl0 13 cm 



triangles: p m — 1.0 x 10 cm 



, squares: 
-3 



Pm = 2.0 x 10 12 cm" 



As expected the agreement becomes better for a larger 
system size (resp. for a smaller macroion density) since 
the theory is constructed formally for vanishing macroion 
density. In addition, we repeated the calculations for 
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p m = 1.17 x 10 13 cto -3 (corresp. to the circles in FigJ|) 
using the PRI method with N = 4 and got the same re- 
sults as in HSG. We now enhance the Coulomb coupling 
by formally reducing the dielectric constant e. For a fixed 
macroion position at Z\ — d m , the force F\ is shown in 
Fig.|| for the parameters of run D. 



= 25 



x x 



X X 



20 40 60 80 100 120 

e 

FIG. 5. Force Fi = F\ ■ e z acting on a single macro-ion 
versus dielectric constant e for a fixed position at Z\ — d m . 
The force is scaled by the (arbitrary) unit Fo = 0.01 . 
The parameters of the system are from run D. The solid line 
is the prediction of PB theory. The crosses are simulation 
results in HSG. 

While the PB theory always predicts a repulsive force, 
the simulation data are in accordance with theory only 
for large e but the force changes its sign for e < 10. Hence 
as expected the theory breaks down for large Coulomb 
coupling where correlation between the counterions be- 
come significant. 

For the same run D, the distance-resolved macroion 
force F\ is shown in Fig.^| for a strongly reduced dielectric 
constant e = 3.9. The simulation data were obtained in 
HSG but confirmed by PRI calculations with N = 4. 

r( 2 ) _ Tp( 2 ) . g z anc [ t ne depletion 



The electrostatic part ' = F} 



part F^ = F^'' ■ e z are shown separately. F z ''' is always 
repulsive and increases with decreasing Z\, at least if 
the macroion is not too close to the surface when the 
counterion depletion between the macroion and the wall 
induced by the finite counterion core is negligible. This 
is an expected behavior, since in general there are more 
counterions close to the walls. The pure electrostatic 

(2) 

contribution, Ff , on the other hand, exhibits a more 
subtle behavior. If the macroion is close to the midplane, 
it is repulsive, then it becomes attractive as the macroion 
is getting closer to the plates. 



?(3) 



(3) 
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but now for run D, e — 3.9, and for 
The crosses are simulation data 



FIG. 6. Same as Fig 
a force unit of Fo 
in HSG for the total force Fi, the squares (resp. circles) are 



simulation data the electrostatic part F^ 2 ' (resp. the depletion 
part F^^.The line is a guide to the eye for the total force. 
The inset shows the effective potential in units of ksT versus 
reduced macro- ion distance Z\/L together with the energy 
barrier AV e //. The solid line is for run D with e — 3.9, 
dashed line is for run E. 

As a function of macroion distance, the total force F\ 
is repulsive, attractive and becomes repulsive again. For 
small separations (which are still larger than the micro- 
scopic counterionic core) the force is dominated by the re- 
pulsive depletion force. Hence the macroion has got three 
equilibrium positions, two of them are stable, namely the 
midplane and a position in the vicinity of the plate. In 
order to extract more information, we have calculated 
the effective wall-particle potential defined by 



Veff(Zl) = - 



2i 



Fx^dh 



(22) 



by integrating our data with respect to the macroion alti- 
tude h. This quantity is shown as an inset in Fig|| One 
first sees that the global minimum is in the vicinity of the 
walls. Furthermore the barrier height A 14// to escape 
from there is about 8/csT. This implies that the time for 
a colloidal particle to escape from the position close to the 
surface is roughly To exp (A14////CBT 1 ) = e 8 r » 3000r 
p8| , pfj[ where tq is a Brownian time scale governing the 
decay of dynamical correlations of the macroion. It can 
also be seen that, for a doubled surface charge (run E), 
the height of barrier increases. 

A similar behaviour occurs for another parameter com- 
binations (run G), see Fig.|^, corresponding to aqueous 
suspensions of micelle-sized macroions. It hence seems 
to be a generic feature of the primitive model for strong 
Coulomb coupling. 
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FIG. 7. Same as FigJ^ but now for run G and for a force 
unit of F = O.lf^. The squares are simulation results 

,14 



for the total force in HSG for a p = 1.19 x 10 14 -^, while 
the circles are for a v — 2.38 x 10 — ^r. The lines are a 
guide to the eye. The inset shows the effective potential in 
units of ksT versus reduced macro-ion distance Z\jL. The 
dashed line is for a v = 1.19 x 10 14 — the solid line is for 
a p = 2.38 x 10 14 -^. 
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FIG. 8. Effective force Fx = Fx ■ e z (circles) and gradient 
of the potential energy Fx — Fx ■ e z (squares) versus reduced 
macro-ion distance Zx/L for run K. The unit of the force is 

z 2 

Fo = 0.1 % . the lines are a guide to the eye. The dashed 
line are data from Ref. [60] 

We note that the barrier height A 14// is about 70fcsT 
which implies a very large escape time. Finally we show 
for a certain parameter combination (run K) which was 
also used in Ref. |30| that the averaged force F\ differs 

from the gradient of the averaged potential energy F%. 
As in Ref. [B0|, the system consists of a single macroion 



in a planar slit of width 5d m , with one charged and one 
neutral wall. Results are given in Fig.||. We conclude 
that the forces behave even qualitatively different. The 
average force Fx is a short-range attractive force which 
becomes repulsive only for touching macroion configura- 
tions. Contrary to that, the force F\ is repulsive up to 
distance about d m /2 from the wall surface. Note that 
our data actually differ from those of Ref. (60| due to the 
early truncation of the Coulomb interaction performed 
there. 



VI. TWO MACRO-IONS BETWEEN PLATES 

We finally consider two equally charged macroions at 
the positions Rx — {Xx,Yx,Zi) and R2 — (X2, ^2,^2) 
confined between plates. A schematic picture is given in 
Fig.|. 




A', 



FIG. 9. Schematic picture for the macroion pair near a 
charged wall of surface charge density a v . For the sake of 
clarity, the position of second wall is omitted. The different 
forces are shown for the case of a mutual attraction. 

In order to reduce the parameter space, we assume for 
simplicity that both macroions have the same altitude 
Z\ = Z2- The distance between the macroion centers is 
R12 =| Rx — R% I where the difference vector Rx2 — Rx — 
R 2 is in the iry-plane. We assume that only one of the 
walls is charged and that the second wall is neutral. This 
gives us the possibility to simulate higher surface charge 
densities. Also for strong coupling, the counterions of the 
two different walls are practically decoupled such that the 
set-up with a single charged wall is not expected to differ 
much from the symmetrical set-up. The total force acting 
on the two macroions can be split into a part pointing in 
z-direction and another contribution pointing along R\2- 
Hence we write Fj = Pf + P± defining = (Pj ■ R 12 



R12/R12 an(J 



F 



Fj-, and Fj 1 



-P! 



for j = 1,2. Clearly, 
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It is instructive to compare these force to the DLVO 
bulk theory which yields 



DLVO 



Z 2 e 2 exp((ri m -R 12 )/Xp) 
eR 12 {l + d m /2X D ) 2 

1 | 1 \ Rl2 

R\2 Xd J Rl2 



(23) 



Here the Debye screening length Xd is given by eq.(18), 
where po corresponds to the counterion number density 
coming only from the macroions, po — ^p m - One can 
also modify the DLVO theory by admitting screening also 
from the counterions stemming from the wall assuming 
they follow a Poisson-Boltzmann density profile. This 
yields the PB force ]|§ 

P[ B = (i? ps )H + {F[ B ) 1 ~ 



(24) 



where we get for the parallel part of the force 
Z 2 e 2 exp(-i? 12 /A D (Z 1 )) 



eR 



1 



+ 



12 
1 



Rl2 



R12 Xd(Zi) 
The perpendicular part of the force is 

Z 2 e 2 A Z3 (Z 1 ) 7o 3 



-jj— ■ (25) 



r>j ,± _ pPB _ 



tan( 7o (l-Zi/L)) 
X cos 2 ( lo {l-Z 1 /L)f z 



(26) 




3 

R I2 /d m 

FIG. 10. Parallel part of the effective force acting onto 
a macroion pair, f\ — f\ ■ R12/R12, versus reduced in- 
terparticle distance Rii/dm- The unit of the force is 

Fo = ^ Z td 2 j x 10 -2 . The parameters of system are from 

run L and for an altitude of macroions of Z\ = 0.6d m . The 
solid line is the bulk DLVO theory, the dashed line is the 
Poisson-Boltzmann theory (^) and the points are simulation 
results in HSG. The statistical error corresponds to the sym- 
bol size. 



The space dependent Debye screening length is 



(27) 



Here p PB (z) and F BB are given by Eqns. ( |l9| ) and (|21 
Contrary to the bulk DLVO force, the PB force has an ad- 
ditional perp endicular part for a pair of particles (second 
term in (pq)). This additional force is attractive. Still 
the total perpendicular force ( p6| ) is always repulsive. 

For the parameters of run L corresponding to weak 
coupling, simulation results for F-j' = Pf ■ R12/R12 are 
shown in Fig.|lO|. 

The solid line corresponds to the bulk DLVO force, 
and the dashed line is the Poisson-Boltzmann result. The 
force is repulsive both in theory and simulation, but the 
theories overestimate the force significantly. As expected 
the Poisson-Boltzmann approach yields better agreement 
than DLVO bulk theory. 

Results for F± for stronger coupling (runs M and N) 
are displayed in Figs.|ll[0. 




FIG. 

Fo = 



11. 



Same as 
x 1(T 3 



1.4 

a 

m 

but 1 



Fig.nDl hut now for run M and 
Simulation results are shown for 



three different surface charges: squares: a p = 0- 



gles: o p = 2.98 x 10 1 



r , circles: a v — 5.95 x 10 



trian- 
r . The 



lines are a guide to the eye. 



For a neutral wall, the interaction force between 
macroions is already attractive. With increasing sur- 
face charge the attraction between macroions becomes 
stronger. Clearly this attraction is neither contained in 
DLVO theory nor in the Poisson-Boltzmann approach 

PD- 
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1.4 1.6 1.8 

2 /d m 

FIG. 12. Same as Fig.hl] but now for run N and for 
Z\ — 0.7 dm. Results are shown for three different surface 
charges: squares: a v — triangles: cr„ 

circles: a p = 2.38 x 10 14 ^. 

y cm* 



, p = 1.19 x 10 14 -^ 

* cm* 



In Figjl^ we fixed the macroion distance and calculated 
f\ and the force perpendicular to the plates, F\ = F^- ■ 
e z versus altitude Z\ for run N. 




1.5 
Z,/d 

FIG. 13. Parallel F\ = Ff ■ R12/R12 ( squares) and 
perpendicular F± — Fj • e z (circles) parts of effec- 
tive force versus reduced altitude Z\/d m for fixed inter- 
particle spacing R12 = 1.2 d m . The unit of the force 

x 10~ 3 , and for the force F\ is 
: . The surface charge density is 



Fl is F = 



x 10 



2.38 x 10 14 - 



The inset shows the effective potential 



in units of fcsT versus reduced macro-ion distance Z\jd m . 

There is attraction. Both the interparticle attraction 
and the wall-particle attraction become stronger in the 
vicinity of the plate. The effective wall-particle interac- 
tion potential for the perpendicular part is shown as an 



inset in Fig.[L3[ Note that the minimum of V e ff is much 
more than twice as large as in the single macroion case 
(compare to inset in Fig]?], solid line). Thus, a pair of 
macroions near a planar surface is more stable than a 
single macroion. This is also evident from the results for 
run N shown in FigJlJ. 



-20 



o -30 



-60 

FIG. 14. 
ing onto a 

* = 

tance R12 / d 








1 



R 12 /d... 



Perpendicular part of the effective force act- 
macroion pair, F\ = F^ ■ e z in units of 

x 1CF 3 versus dimensionless interparticle dis- 

The parameters of system are from run N and 
the altitude of macroions is fixed to Zi = 0.7 d m . Simulation 
results are shown for two different surface charges: squares: 



1.19 x 10 1 



circles: a v = 2.38 x 10 1 



Again there is attraction towards the plate for varied 
R12 and fixed Z±. The attraction becomes stronger if 
the interparticle distance is decreasing. This shows that 
the attraction between the wall and a single macroion 
discussed in chapter V is stable and even enhanced if 
more macroions are close to the wall. This leads us to the 
final conclusion that the macroions will assemble on top 
of the surface forming two-dimensional colloidal layers. 



VII. CONCLUSIONS 

We have simulated the effective force between 
macroions confined in a slit-geometry. An effective at- 
traction was found for strong Coulomb coupling. In par- 
ticular, the effective potential of a single macroion con- 
fined between two parallel charged plates was found to 
have two stable minima where the total force vanishes: 
the first is in the mid-plane, the second close to the walls. 
This result was confirmed for two macroions. In this case 
the attraction towards the walls was even stronger than 
for a single macroion. Our most important conclusion 
is that the attractive force will result in two-dimensional 
colloidal layers on top of the plates. As the depth of 
the attractive potential is larger than fcgT, these layers 
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possess a large life-time with respect to thermal fluctu- 
ations. The layers should be crystalline as the interpar- 
ticle interaction is also attractive. This can explain at 
least qualitatively the long-lived metastable crystalline 
layers found in recent experiments on confined samples 
of charge colloidal suspensions j3^J|^] . 

We want to add some remarks: First, our parameters 
are actually different form those describing the experi- 
ments. The main difference is the high surface charge 
of the glass plates within an area spanned by a typical 
macroion separation distance. Such a system cannot be 
simulated since it requires a huge number of counterions 
in the simulation box. We have mimicked the high sur- 
face charge by dealing with a small dielectric constant, 
but strictly speaking this corresponds to a different sys- 
tem. Second, the mechanism of our attraction is similar 
to that proposed recently by us in the bulk case f25fl , ft 
only occurs for strong coupling with divalent counterions 
and is short-ranged. In this respect, it behaves differ- 
ent than in experiment where the attraction was long- 
ranged. We emphazise again that the depletion force is 
crucial in the strong coupling parameter regime. Third, 
our computer simulation data were tested against simple 
DLVO- or Poisson-Boltzmann-type theories. It would be 
interesting to use them as benchmark data for more so- 
phisticated theoretical approaches which predict attrac- 
tion as e.g. the density functional perturbation theory re- 
cently proposed by Goulding and Hansen [^|. Finally, in 
our simulations, we neglected any impurity or added salt 
ions. Their inclusion increases substantially the number 
of microscopic particles and would lead to more extensive 
simulation. A further challenge would be to incorporate 
image charges properly into the model which requires a 
non-trivial extension of our approach. 
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APPENDIX A: DEFINITION OF FORCES IN 
HYPER SPHERE GEOMETRY 

We shortly present here some technical details in HSG 
within the primitive model. For more details, we refer to 
Ref. || . 

The charged hard spheres are confined on the surface S3 
of a 4D hypersphere. Without confinement, in the bulk, 
the whole surface of the hypersphere is accessible to the 
particles. Since S3 is compact, the total charge in a closed 
space must be equal to zero. We define a pseudocharge as 



the association of a point charge q\ located at the point 
say M, and a neutralizing background of total charge 
—q. L . The position of pseudocharge is specified by the 4D 
spherical coordinates (a, 6, ip) (see Fig. fL5|) . 



North pole 




South pole 



FIG. 15. Schematic view of the hypersphere (projected to 
three dimensions) illustrating the angular coordinate a. 

Then the Cartesian components of the unit vector 
u(M) = OM/R reads 

U\ = sin a sin 9 cos ip, u-i = sin a sin 9 sin ip, 

U3 — sin a cos 9 , u^ — cosa (Al) 

Here R is the hypersphere radius, O is the center of the 
hypersphere and the angle a determines the distance Ra 
from north pole, see Fig. |l5|. The distance between two 
pseudocharges (at point M) and qj (at point N) is 
measured along the geodesic line joining these points 

r M N = Rl, (A2) 

where 7 is the angle between vectors OM and ON, 

OM ■ ON 



7 = arccos ■ 



R 2 



(A3) 



The Coulomb force F^j between pseudocharge and 
pseudocharge qj is 



tR 2 



F, , =-^fcot 7 +^)e 7 (Ar). 



. 2 

sin 7 



(A4) 



Here e 7 (N) denotes the unit vector tangent to geodesic 
MN at point N 



1 



sin 7 



-u(M) + cot 7 u(N) 



(A5) 



For short distances tmn, there is the hard core repulsion. 
A hard sphere of diameter di (i = c, to) centered around 
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the point M on S3 is denned as the set of points Mq 



such that RjMM 



R arccos QM-OMn 



< di/2. Thus, 



the hard sphere potential between two pseudocharge is 
defined by 



Ui 



00 




. 2R ■ 

otherwise. 



if 7 



hJ 



(A6) 



Let us now consider the mixed hard sphere system to 
be confined between two charged walls. On S3 geometry 
it corresponds to the two charged lamellae, parallel to 
each other, localized symmetrically on opposite side of 
the equator (see Fig.(p^)). 



North pole 
lamellar ] 




South pole 



FIG. 16. Schematic view of the hypersphere (projected to 
three dimensions) illustrating a situation with two parallel 
charged interfaces. 

They result from a conical section of the hypersphere, 
with angular apertures equal to ajv for north lamellar 
and as — ft — cvjv for south lamellar. The area of each 
lamellar is S p — 4nR 2 sin 2 a jy and the volume confined 
between lamellae is given by 



V(cxn) — irR 3 (2aiy ~ sin2aAr) 



(A7) 



Then the separation distance between lamellar is D = 
R(t: — 2a n)- For the symmetrical case considered in this 
paper, when both lamellar are charged equally with sur- 
face charge density a p , the charge electroneutrality of 
simulation cell together with eq. (Qj) requires 



2Q P = 0. 



(A8) 



Here Q p = a p S p is the net charge of one lamellar, N p is 



the number of counter- ions coming from planes. 

Finally, we give the definition of interaction forces. The 
lamellar-lamellar repulsion force is 



Ql 

irR 2 



cot a jv 



sin ajv 



(A9) 



The ion-ion repulsion (i = j) and attraction (i ^ j) force 
is 



7T — 7 

~ 2 

sin 7 



(A10) 



here 7 is given by (A3) and the particles i and j are at 
points M and N, i,j = m,c. The ion-lamellar repulsion 
(i = m) and attraction (i = c) force is 



F„ 



Qp Qi 



enR 2 

for the upper lamella and 

Q P q, 



cot cvi 



Fvp ~ enR 2 



cot as 



sin 2 a* 



sin 2 a, 



(AH) 



(A12) 



for the second lamella. The direction of the forces is 
always along the geodesic line. 

The lamellar-ion hard core repulsion becomes 



00 for 04 < aN + ck/2R, i = c,m 
otherwise 



(A13) 



APPENDIX B: EQUATION OF MOTION FOR 
SINGLE COUNTER-ION IN HYPERSPHERICAL 
GEOMETRY 

In this section we translate Newton equation of motion 
onto HSG for a counterion in an external electrical field 
created by charged planes, other counterions and fixed 
macroions. First of all, let us define the displacement of 
counterion at point M on S3. The differential dOM of 
vector OM (M remaining on the surface of the hyper- 
sphere) is 

d OM = Rdae a + R sin adO eg + R sin a sin 6 e v (Bl) 

where the unit vectors [e a ,e$,e v ] given by 

e„ = (cos a sin 9 cos <p, cos a sin 9 sin tp, 
cos a cos 9, — sin 9), 



(cos 6 cos <p, cos 6 sin ip, — sin 0,0), 



(— sin (p, cos ip, 0, 0). 



(B2) 



constitute an orthogonal frame at point M. For the ki- 
netic energy in terms of the variables (a, 9, tp) we get 



T 



mv m 



( dOMX 
~2~ ~ ~2 I dt J 



in 



(i? 2 d 2 + R 2 sin a 2 9 2 + R 2 sin a 2 sin 9 2 ip 2 ^j (B3) 
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For the potential energy we have relations 

dU dU 
- — = RF a , — — = RsmaFg, 
oa oO 

dU 

— — — — R sin a sin 9 F a 

dip 



(B4) 



Now let us put Eqs.QB3|) and ( |B4| ) into Lagrange equa- 
tions 



N r . 



d dT 8T _ dUij 
dt dvi dxi i— 1 



3=1 



dxi 



E 

k=l 



ki 



dxi 



E 

1=1 



dUu 
dxi 



(i = a,6,(p and Xi = i) (B5) 



where on the right hand side the first term arises from all 
counterions, the second term is the macroionic attraction, 
the last term is due to plane attraction. The Newton 
equations of counterion motion reads as follows 



a = sin a cos a ( 9 Z + sin 2 9ip 2 



1 



m R 



-(sm*a0) 
dt\ J 



•2 • • 2 8111 a 

= sin a sin 9 cosd <p H Fg 

m R 



( • 2 -2 

— sm a sm 
dt v 



sin a sin# 
m R 



(B6) 



(B7) 



(B8) 



We note that (F a , Fg, F v ) are the components of to- 
tal force arising from all other counter-io ns, p lane s and 
macro-ions. The equations of motion Eq. (|Bg )-(Bj ) have 
been solved in a way similar to that described in ]3l| . 

The spherical components (F a , Fg, F^) of vector F are 
connected with the 4D force eq.( A4) by the following 
relations 

F a = Fi cos a sin 9 cos tp + F2 cos a sin 9 sin ip 

+ F 3 cos a cos 6 — Fa sin 9, 

Fg = Fi cos 9 cos ip + F2 cos 9 sin ip — F3 sin 0, (B9) 

F v = —Fi sin + F 2 cos <p. 
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